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ABSTRACT 

Numerical  methods  are  applied  to  one-dimensional  unsteady 
reaction-diffusion  equations  to  seek  propagating  wave  solutions. 
These  equations  describe  flame  propagation  in  certain  combus- 
tion systems  if  constant  pressure  combustion  is  assumed,  with 
Lewis  number  =  1  and  if  a  Lagrangian  coordinate  transformation 
is  introduced.   It  is  shown  that  the  steady  flame  speed  is 
invariant  under  this  coordinate  transformation.   Model  reaction- 
diffusion  equations  which  admit  traveling  waves  as  exact  solutions 
are  formulated  and  one  of  these  is  solved  with  twelve  different 
numerical  integration  schemes  in  a  test  problem.   The  diffusion 

terms  are  differenced  using  the  explicit  methods  introduced  by 

3    2 
Saul'yev,  which  are  shown  can  be  formally  accurate  to  0(At  /Ax  ). 

Both  implicit  and  explicit  techniques  for  the  reaction  rate 
terms  are  tested.   The  results  of  a  study  with  these  schemes 
indicates  that:   numerical  diffusion  can  reduce  accuracy  signifi- 
cantly, that  numerical  dispersion  truncation  errors  can  reduce 
accuracy  if  they  are  sufficiently  larae ,  and  that  an  accurate 
representation  of  the  reaction  rate  in  the  difference  equations 
is  important  to  retain  overall  accuracy.   In  addition  to  the 
above  test  problem,  one  of  the  explicit  methods  is  applied  to 
an  ozone  decomposition  flame  computation.   An  adaptive  grid 
procedure  is  also  implemented.   The  results  show  good  agreement 
with  the  fourth  order  accurate  results  of  Margolis  which  indicates 
that  a  more  efficient  lov/er  order  numerical  method  can  be 
sufficiently  accurate  for  practical  computations. 
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1.  INTRODUCTION 

Detailed  mathematical  models  are  now  being  applied  to  prac- 
tical problems  in  combustion.   Attention  is  being  given  to  models 
which  include  comprehensive  elementary  chemical  kinetic  and  tur- 
bulence mechanisms.   This  trend  has  increased  the  number  of 
equations  to  be  solved  and  has  stimulated  research  in  numerical 
methods  for  the  conservation  equations  of  reactina  gases. 

In  this  paper  we  discuss  results  for  one-dimensional  uncon- 
fined  deflagration  wave  propagation.   The  problem  has  been  for- 
mulated in  such  a  way  that  the  momentum  and  overall  mass  conser- 
vation equations  need  not  be  considered.   This  simplification 
focuses  attention  on  the  treatment  of  reaction  rate  and  diffusion 
terms  in  the  governing  equations. 

The  presence  of  reaction  rate  terms  can  cause  timestep 
restrictions  in  computations  of  combustion.   Each  chemical  species 
evolves  with  its  own  characteristic  time  scale,  and  the  numerical 
timestep  may  be  required  to  be  comparable  to  the  smallest  of 
these  time  scales  for  reasons  of  stability  and/or  accuracy. 
This  problem  has  led  many  researchers  to  use  implicit  solution 
techniques  in  combustion  calculations. 

Reaction  diffusion  equations  have  also  been  solved  following 
this  approach.   This  includes  methods  in  which  a  system  of  coupled 
ordinary  differential  equations  (O.D.E.)  are  solved  using  stiff 
O.D.E.  routines  (method  of  lines,  Bledjian  [1],  Margolis  [2] ; 
operator  splitting,  Otey  and  Dwyer  [3]);  and  methods  in  which 


block  tridiagonal  matrices,  which  result  from  an  implicit  treat- 
ment of  the  diffusion  terms,  are  computed  (Lund  [4]) . 

In  this  work  we  formulate  and  apply  an  explicit  method  to 
a  one-dimensional  ozone  decomposition  flame  calculation.   The 
motivation  for  the  use  of  an  explicit  method  is  the  potential 
for  increases  in  efficiency  and  accuracy.   These  benefits  are 
important  considerations  in  flame  propagation  problems  which 
involve  a  large  number  of  participating  chemical  species.   The 
numerical  method  which  we  used  for  the  ozone  flame  problem  does 
not  have  the  diffusional  stability  restriction  but  the  numerical 
timestep  is  limited  by  the  Lipschitz  timestep  constraint, 
which  arises  from  the  reaction  rate  terms.   However,  the  results 
of  Section  4  of  the  paper  show  that  the  timestep  of  this  explicit 
method  is  nevertheless  comparable  to  the  timestep  used  in  the 
implicit  method  of  Lund  [4] . 

We  also  compare  the  results  of  several  other  numerical 
schemes  which  use  implicit  methods  for  the  evaluation  of  the 
reaction  rate.   In  that  study  however,  we  employed  a  single 
non-linear  reaction-diffusion  equation  (which  admits  exact 
solutions)  as  a  model  problem. 


2.  GOVERNING  EQUATIONS 

The  general  set  of  equations  which  describe  one-dimensional 
flame  propagation  consists  of  the  conservation  and  state  equa- 
tions for  reactive  flows  which  can  be  written  as  follows 
(Williams  [5]). 
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The  first  four  equations  are  the  overall  mass,  momentum, 
energy  and  the   N   species  conservation  equations  for  the   N 


chemical  species  having  mass  fractions 
mass  averaged  gas  velocity. 


(k) 


is  the 


pressure,  temperature  and  density.   In  Eq.  (5),  the  thermal 

equation  of  state,  it  has  been  assumed  that  the  mixture  behaves 

(k) 
as  a  perfect  gas.   Here  W     is  the  species  molecular  weight. 

Fick's  law  has  been  employed  to  describe  species  diffusion 

effects  with  species  diffusion  coefficients  ^  .      to      is 

the  mass  rate  of  creation  or  depletion  of  species   (k)  ,   and 

h^  '   is  the  heat  of  formation  of  species   (k)   at  temperature 


Thermal  radiation,  body  forces,  Soret,  Dufour,  and  pressure 
gradient  diffusion,  and  bulk  viscosity  effects  have  been  assumed 
to  be  negligible. 

Model  Formulation 

In  unconfined  flame  propagation  the  combustion  wave  separates 
the  unburned  gases  at   x  =  °°   and  the  burned  gases  at   x  =  -°° 

as  depicted  in  Figure  1.   The  unburned  mixture  has  constant  density   p 

(k) 
temperature   T    and  species  mass  fractions   Y     and  velocity  u   =0. 

In  the  burned  gases  p  =  p^,  T  =  T   and  Y^^'''  =  ^V  •  However,  the  quan- 
tities at   X  =  -«>   and  the  wave  speed,   S  ,   are  unknown  in  general 
and  are  to  be  computed. 

The  governing  equations  can  be  reduced  to  a  set  of  reaction- 
diffusion  equations  if  simplifying  assumptions  are  introduced  and 
if  a   coordinate  transformation  is  made.   The  constant  pressure 
assumption,  P3^=Pq  valid  when  the  flame  velocity  is  very  much  less  than 
the  sound  speed,  allows  the  momentum  equation  (2)  to  be  drooped. 


Another  consistent  approximation  is  to  neqlect  the  terms  involving  the 

2 

kinetic  energy  ('^^u  /2)  and  its  derivatives  in  Eqs .  (3)  and  (6)  (see 

Williams  [5]),   The  convective  terms  in  (3)  and  (4)  can  be  eliminated 

Lagrangian 
by  a  /coordinate  transformation  from  the  (x,t)  plane  to  the  (x,t)  plane 
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and  satisfies  Eq .  (1)  identically.   Equation  (9)  may  be  rearranged 

using  (8)  and  (6)  (after  summation  over  (k)  and  multiplication  by 
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h    )  to  give 
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(k) 
Finally,  with  the  further  assumptions  that  ^or   all  k,   c     =  c   , 

^(k)  ^^    where   c    and  ^  are  constants,  and  that  p^  =    ^/c    (Lewis 
number  =  1)  we  have  the  reaction  diffusion  equation  system 
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Supplemented  by  the  equation  of  state 

pTZY^'^Vw^'^^   =   p/.^   =   const. 


(12) 


and  the  boundary  conditions  (which  enforce  zero  diffusion  of  heat 
or  mass  at  ±°o)  . 


9Y 


(k) 


9x 


9t 
9x 


9Y 


(k) 


=   0 


9x 


(13) 


for   k  =  1 , . . . , N  . 

Conventionally,  the  gas  velocity,   u  ,   and  the  wave  speed 
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S,  is  determined  by  integration  of  Eq .  (7a, b),  once  Eqs.  (10) 

to  (13)  have  been  solved.   However,  notice  that  in  steady  propagation 

the  wave  speed  is  invariant  under  the  transformations  (7a, b) 

This  is  because  in  steady  wave  propagation 

|-+s|-  =   ^+s—  =   0   ,  (14) 

3t      3x     g:^      g- 

and  therefore  (4)  and  (8)  may  be  integrated  (subject  to  Eqs. 
(13)  which  are  the  same  in  both  coordinate  systems)  to  show 
that 


p(u  -  S)Y^^' 


(15: 


(k)   _ 
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Now  equation  (1)  shows  that   p (u  -  S)  =  const. p^S  ,   and 

the  use  of  (7a)  in  (15)  and  a  comparison  with  (16)  verifies 
that   S  =  S  .   This  observation  simplifies  the  determination  of 
the  steady  state  flame  speed  from  computational  results  obtained 
in  the   x,t   plane.   In  summary  then,  Eqs.  (10),  (11)  and  (13) 
are  written  in  the  generalized  form 
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8  f^(k)  du^^^ 


9t       8x [        3x 
subject  to  the  boundary  conditions 

We  now  consider  the  case   K  =  1   and   d^"^^  =  d  =  constant. 

u  =  u      could  be  taken  to  represent  non-dimensional  temperature 
such  that   0  ^  u  ^  a  ,   where   a   is  the  adiabatic  flame 
temperature.   This  simplification  of  Eqs .  (10)  and  (11)  can  be 
shown  to  be  valid  for  the  case  of  a  global  irreversible  one-step 
chemical  kinetic  mechanism  with  a  stoichiometrically  premixed 
mixture  of  reactants. 

In  this  case   F(u)  (>^0)   is  usually  given  by  a  non-linear 
Arrhenius  expression  (see  Eq .  (40)).   Spalding  [6]  proposed 
a  simpler  algebraic  expression 

F(u)  =  3u'^(a  -  u)""  ,  (IG) 

which  approximates  the  Arrhenius  function  for  large  positive,  n  , 
but  no  exact  solutions  have  previously  been  found  for  this  model. 

It  is  noted  here  that  exact  travelling  wave  solutions  can 
be  found  for  n=2  and  m=l  or,  more  generally,  for  F  given  by 

F(u)  =  6u''^^(a  -  u"")  ,    m  >  0   .  (19) 


The  solutions  are  given  by 

mS ,   _ . .   .  , 
r  r  -5— tx-St)^-|l/m 

u(x  -  St)   =    a/  1  +  e"^  ,  (20) 


and  they  correspond  to  the  class  of  stable  minimal  wave  speed  solu- 
tions, with  exponential  decay  for  x  -^  «>  ,  which  are  discussed  by 
Lin  [16]  . 
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In  Eq.  (20)  the  wave  speed  is 

S   =   a(3d/m+l)^ 

An  additional  model  equation  is   n  =  m  =  1   in  Eq .  (18) . 
This  problem  (Fisher's  Equation)  has  been  analysed  theoreti- 
cally by  Kolmogoroff  et  al.  [7]  and  numerically  by  Gazdag  and 
Canosa  [8]  in  studies  which  show  the  existence  of  a  minimum 
wave  speed.   McKean  [9]  has  shown  that  the  final  steady  state 
wave  speed  is  determined  by  the  nature  of  the  initial  data. 

In  particular,  for  initial  data  which  decays  exponentially  ahead 
of  the  wave,  i.e. 

u(x,0)   =   ae"'^^      (x  ^  -)  (21a) 

where   a   and   b   are  constants,  the  asymptotic  steady  state 
wave  speed  is 

S   =   b  +  1/b     .  (21b) 

Ablowitz  and  Zeppetella  [10]  have  found,  for  the  particular  wave 
speed  S  =  5//6"  ,  that  an  exact  travelling  wave  solution  of  the 
Fisher  Equation  is  given  by  (here   a  =  B  =  d  =  1) 

u(x  -  5//6  t)   ^   [l  +  e(^-^/^^)/^^)  "'    .  (22) 
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3.  NUMERICAL  METHOD 

We  describe  an  explicit  numerical  method  for  reaction- 
diffusion  equations  which  will  be  shown  to  be  accurate  to 

3   2 
0(At  /h  )   under  certain  conditions.   The  method  is  given  in 

Eq.  (23)  below.   Variations  of  the  method,  and  particularly 

several  implicit  treatments  of  the  reaction  rate  terms,  are 

obtained  by  varying  the  parameters  y,    5,    a,    q,    ^      in  the  equations. 

Discussion  of  the  implicit  methods,  6  7^  0,  is  however  postponed 

until  the  next  section.   In  this  section  we  also  describe  an 

adaptive  grid  method  which  is  used  in  the  ozone  flame  computations 

to  be  discussed  in  section  4. 

We  consider  a  finite  difference  scheme  in  which  the  diffu- 
sion terms  in  Eq,  (17)  are  differenced  using  the  explicit  method 

first  introduced  by  Saul'yev  [11].   He  was  only  able  to  show  that 

2    2 
the  numerical  method  was  accurate  to   0(At  /Ax  ).   The  finite 

difference  approximation  to  the  continuum  solution   u"  .  'v  u    (x  .  ,nAt; 

'^  i3    —  D 

is  computed  at  grid  points   x .   from  the  set  of  equations 


J k^   =    _1 u     Ln  ^^n  Ln  _^n 


At  -    h^h^h     ihl[^-k,j  +  l-'^k,jj["k,j  +  l""k,j 


-  ^2[A<]<U]  (vrVj-i)-(i-Y)  (^k,i^^k,:-i)  K-,r"K, j-i)' 


n-l- 


6F^(V.)    +    aF^(u';^j^g^^)     +    BF^iu,^.]  • 


(23a) 
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1   k,i 


At 


E^{-sK,i-^^.i-i)K,.-"i:,i-i) 

+  6F^(W,)  +  aFju5_,^q^,)  +  PF^K;^)  (23b) 


u"*i  =  (V.  +  WA/2  (23c) 


where 

^1  ^  ^j  -  ^j-1  '  ^2  =  ^j  +  1  -  ^j  '•   h  =  h^  +  h2 

^1  ^  ^i  -  ^i-1  '•  ^2  ^  ^i+1  -  ^i  '  h-  =  h!^  +  h^ 

j=l,2,...,M;  i=M-j+l 

k=l,2,...,K;  Jl  =  l,2,...,K 

an(i   M   is  the  number  of  computational  points.   In  the  present  dis- 
cussion we  take 

6  =  q  =  0 

a  -  3/2  ;   6  =  -1/2  (23d) 

In  this  case',  Eq .  (23a, b)  consist  of  two  predictor  sweeps  of  the 
mesh  (in  opposite  directions)  followed  by  a  corrector  step  in 
Eq.  (23c) ,  with  the  reaction  terms  being  given  by  an  explicit 
extrapolation  formula.   Computationally,  the  operation  count  is 
similar  to  that  of  the  back  substitution  steps  in  the  LU  decom- 
position of  the  matrices  in  certain  implicit  methods. 
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The  difference  equations  (23)  are  still  implicit  in  the 
diffusion  coefficients  d?,     .  .   This  difficulty  is  avoided  if 

it  is  assumed  that  d^    .    '^   d,      .      or  if   d,     can  be  computed 

from  quantities  whose  values  are  already  known  at  time  level 
n+1  .   We  will  assume  the  latter  to  be  the  case.   In  addition, 
Eqs.  (23)  are  also  implicit  in  the  boundary  conditions, 
but  this  is  not  of  concern  in  wave  propaaation  problems  due  to 
the  specification  of  Neum.ann  boundary  conditions  as  in  Eqs.  (13) 
(Roache  [12] )  . 

We  now  show  that  the  numerical  method  represented  by  Eqs.  (2  3) 
has  leading  terms,  given  by  a  formal  truncation  error  analysis,  which 
are  similar  to  those  of  the  second-order  accurate  Crank-Nicolson  [13] 
method.   Notice  that  the  Crank-Nicolson  m.ethod  is  obtained  by  replac- 
ing V.  and  W.  by  uf^"*""!"  in  Eqs.  (23)  with   Y  =  5  =  a  =  1/2   and 

q  =  6  =  0  .   In  the  present  example  we  use  an  equally  spaced  mesh 

^1  "  ^2  "  ^/^  • 

Expanding  terms  in  (23)  in  a  Taylor  series  about 

the  continuum  solution   u    (x.,  (n+1)  At)   shows  that  the  continuiom 
analog  of  the  difference  equations  is 

At 


F(u)  =  ^[u^  -  (du^)^  -  F(u))^ 
fl->)(d"x)xt-  H^(^-^)x 


At 
2 


+  0[Yh(V-W)    ,At  ,h^,Ath)  (24) 
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where  the  superscript   (k)   has  been  dropped  for  simplicity  and 
the  subscripts   t   and   x   indicate  differentiation.   The  second 

and  third  terms  on  the  right  hand  side  of  Eq.  (24)  and  the 

2 
0(h(V-W)    )   term  (which  will  be  seen  to  be   0(At  )  ),  represent 

a  departure  from  the  first  order  Crank-Nicolson  truncation  error 
form.  The  term  involving  (V-W)  may  be  estimated  by  subtract- 
ing the  Taylor  expansion  of  Eq .  (23b)  from  that  of  (23a)  giving 


with  this  result,  a  regular  perturbation  analysis  with  small  para- 
meter At  shows  that 

V  -  W  ^  -Mi!X(du^)^  +  ^[(au^)^^    -    Y[d(du^)^JJ  t   H.O.T. 

Substituting  this  equation  in  Eq.  (24)  gives 

"t  -  (^^x^x  -  ^(^)  =  I^K  -  (^^x^x  -  F(u)]t 

+  |t(2Y2r  -  (1  -Y))  (du^)^^ 
2  2, 


At  Y  r 


P^^-x^ttx  -  ^^^(^^Vxt^xx^  ^«-^-^-   ( 
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where   r  =  dAt/(h/2)   .   From  this  result  it  is  seen  that  if  d  =  const 

and  

/l+8r  -  1 

Y   =   (26) 

4r 

3   2 
the  scheme  is  accurate  to   0(At  /h  )  .   Formal  second  order 
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accuracy  only  occurs,  under  the  additional  constraints 

F(u)  =0   and   Y  =  1   which,  in  view  of  (26)  implies  that   r  ~  0  . 

The  condition  for  diffusional  stability  of  the  method  in 
Eqs.  (23)  is  (Saul'yev  [11]) 

p   >   2(1-Y)  (27) 

and  it  is  thus  unconditionally  stable  for   y  >  1  .   However, 
the  Lipschitz  time  step  constraint 


At   <    ^°"^^ 


I  9F,  (U„  .) 
max     k   £,] 

'     m,] 


(2i 


is  present,  and  recently  Hoff  [14]  has  shown  that  this  can  make 
an  unconditionally  stable  scheme  for  the  diffusion  operator  only 
conditionally  stable  for  the  combined  reaction-diffusion  system. 

Finally,  in  this  section  we  describe  the  variable  timestep 
and  adaptive  grid  procedures  which  were  used  in  the  ozone  decom- 
position flame  computations.   The  allowable  timestep  was  found, 
not  from  Eq .  (28) ,  but  rather  by  requiring  that  the  species  mass 
fractions  remain  positive  at  every  computational  point.   If  negative 
mass  fractions  were  found,  the  timestep  was  halved  and  that  compu- 
tational cycle  was  repeated;  after  20  such  computational  steps  the 
timestep  was  doubled.   A  three  time  level  explicit  scheme  (Vlllb  in 
Table  I)  was  used  in  those  computations.   Interpolation  (or  extra- 
polation) between  time  levels,  with  a  change  in  timestep,  was 
avoided  by  adjusting  the  parameters   a   and   g   in  Eqs.  (23a, b). 
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The  parameters  were  changed  such  that  the  first  order  truncation 
error  of  the  reaction  term  remained  equal  to   -y-  F   (see  Table 
1;  scheme  Vlllb)  where   At   now  corresponds  to  the  new  timestep. 

The  adaptive  grid  method  is  similar  to  that  used  by  Lund  [4] 
but,  because  the  mesh  function  is  prescribed  analytically,  it  is 
simpler  to  implement.   The  non-uniform  mesh  was  aenerated  using 
parameters  from  the  solution  of  an  associated  eigenvalue  problem. 
A  new  mesh  was  generated  when  changes  in  the  numerical  solution 
occurred  beyond  a  prescribed  tolerance,  and  the  solution  variables 
were  interpolated  onto  the  updated  mesh.   The  details  are  as 
follows . 

If   z  (x)   is  a  continuous  variable  which  represents  the 
mesh  spacing  at  point   x  ,   a  variational  problem  is  formulated 
which  minimizes  the  integral 


i<->    '    l"i^fall'    '  <"> 


Subject  to  the  constraint 


„.   .   I 


dx 
z 


(3o: 


Minimizing  the  integral  in  Eq.  (29)  minimizes  the  change  in  mesh 
spacing  across  the  domain  where  x,  and  x..  are  fixed  boundary 
points,  while  Eq.  (30)  conserves  the  number  of  zones   M'  =  M  -  1  . 
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Equations  (29)  and  (30)  lead  to  the  Euler  equation 


(31) 


which  is  to  be  solved  subject  to  the  boundary  conditions 


(32) 


Equation  (31)  may  be  readily  integrated  to  give  the  parabola 


^U.c,,^.§- 


where  the  constants  of  integration  c,  and  c^   may  be  determined 
from  (32) .   However,  since  the  desired  mesh  spacing  at  the  boun- 
daries (z,  and  z„)  is  not  known,  c,  and  C2  were  found  using  an 
alternate  procedure.   The  position  where  the  largest  gradient  of 
9t' 


temperature 


8x 


occurs  is  denoted  by  x  .  .   By  reauiring  that 
-"   min 


the  value  of  the  mesh  function,  z  .  ,  coincides  with  this  point, 

min 


the  constants  c,  and  c^  are  determined,  i.e. 


z(x)  =  -7- X  -  x^.     +  z^. 

4z  .        mm      min 


(33) 


where 


AT 


(34) 
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c  is  a  specified  constant  which  is  related  to  the  number  of 
zones  in  the  region  of  large  temperature  gradient,  and   AT   is 
the  difference  between  the  maximum  (T2)  and  minimum  (T-,)  values 
of  temperature  in  the  domain. 

The  eigenvalue   A   was  determined  by  iteration  from  the 
equation 

_^   x=x 

2 


M'/A 
2 


M 
x=x. 


35) 


Eq.  35  is  obtained  by  integrating  Eq .  (30)  with   z (x)  given 

in  (33)  .   The  special  case   A  =  0   corresponds  to  an  equally 

2 
spaced  mesh,  while   A  -    (27t/M')    corresponds  to  the  case 

x^  ^  °°  ,   X,  ^  -°°  .   Equation  (35),  with  the  upper  limit   x 
replaced  by   x.   and   M'   replaced  by   j-1   then  served  to  define 

the  new  mesh,  for   j  =  2,...,M-l  .  The  mesh  was  redefined  if  the 

new   x  .    and   z  .    violated  either  of  the  two  inequalities 
mm        mm 

o 
x.-z.    <x.    <x.+z. 

min    mm   —   mm   —   min    min 

o 

^min   ^         ^   „  o 

—^ <   z  .    <   2z  . 

2     —   mm   —    mm 

where   x  .  ,z  .    defined  the  mesh  at  the  previous  remap, 
mm   mm 

The  solution  variables  (e.g.  the  temperature   T^   at   x^) 

were  mapped  onto  the  new  mesh   x.   with  the  non-linear  interpolation 
formula 
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where 


T,+h .T„ 
1  3  2 
1+h. 


(36a) 


1-a  a 


'i+1 


and   a 


for   x^  £  X .  <^  ^i  +  i   ^"^   T^  7^  '^l''^2 


This 


interpolation  is  exact  for  the  model  Eq .  (20)  for  the  case 


1  .   For   T 


T„  (£=^1,2)   we  used  instead 


i+1 


i+1 


1-a 


(36b) 


which  is  exact  for  exponential  functions.   Equations  (36)  were 
chosen  in  anticipation  of  solutions  which  have  aDnroximately 
exponential  decay  outside  of  a  narrow  flame  zone. 

Finally,  mass  conservation  was  enforced  on  the  new  mesh  by 
adjusting  the  N  species  mass  fractions  at  each  grid  point  x- 
such  that 

Y(k) 
Y(k)   ^     1 
J 


N 

I 
k=l 


is  the  SDecies  mass  fraction  obtained  from  the 

J 

interpolation  Eqs .  (36).   In  practice,  the  correction  was 
found  to  be  of  the  order  of  the  machine  round-off  error  after 
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the  first  few  mesh  remap  cycles,  although  the  correction  was 
Significant  for  the  discretized   Y.     on  the  first  specified 
mesh  at   t  =  0  . 
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4.  DISCUSSION  OF  RESULTS 

We  first  present  the  results  of  a  parameter  study  of 
the  numerical  method  described  in  Eqs.  (23  a,  b,  and  c).  The 
parameters  used  are  listed  in  Table  I.   The  various  schemes 
include  implicit  and  explicit  methods  for  the  evaluation 
of  reaction  rate  terms.   A  single  non-linear  reaction- 
diffusion  model  equation  was  used  in  these  studies.   We  show 
later  our  computational  results  for  a  coupled  system  of 
reaction  diffusion-equations  describing  the  propagating 
ozone  decomposition  flame. 

The  Fisher  Equation,  chosen  as  a  model  equation,  can  be 
written  in  non-dimensional  form  as 

u^   =   du   +  u(l  -  u)  (37; 


which  yields  travelling  wave  solutions  with   0  5_  u  <_  1  . 
Notice  that  the  quadratic  form  of  the  non-linear  reaction  rate 
term  simplifies  the  computation  of  implicit  methods   ( 6  7^  0 
in  Eqs.  (23a, b))  since  they  can  be  solved  explicitly. 

We  set   d  =  1   and  used  an  equally  spaced  mesh  with  the 
domain   -50  ;f_  x  <_  400.    This  domain  was  found  to  be  suffi- 
cient to  eliminate  boundary  effects.   The  initial  condition 
U(x,0)   was  given  by  the  exact  solution  Eq .  (22),  which  is  a 
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travelling  wave  with   S  =  5//6"  and   u(°°,t)  -    0  ,   u(-°°,t)  =  1  . 

In  the  results  shown  in  Fig.  2,  however,  this  initial 
condition  was  modified  slightly.   This  figure  shows  the  numeri- 
cal solution   U(x)  ,   for  various  times  up  to  time   t  =  20 
using  the  standard  explicit  scheme  VI  in  Table  I  on  a  fine  mesh 
with   Ax  =  0.25  ,   At  =  0.01  .   The  initial  data  was  given  by 
Eq.  (22)  for   U  >  5xl0~^  ,   but  for   U  <_   5xl0~^   was  given  by 
Eq.  (21a).   The  constant   b   in  (21a)  was  set  equal  to  0.127, 
i.e.   b  +  1/b  2i  8  ,   and  the  constant   a   was  fixed  by  matching 
the  initial  data  at  the  change-over  point.   Notice  that  with 
these  conditions  Eq .  (21b)  would  predict  a  change  in  the  wave 
speed  from  its  initial  value  of   5//6   to  a  final  speed   S  ^i  8  . 

The  results  shown  in  Fig.  2  indicate  that  the  wave 
propagates  to  the  right  with  its  structure  preserved  up  to 


the  wave  appears  to  settle  to  a  new  structure.   The  corresponding 
change  in  the  wave  speed  is  shown  in  curve  1  of  Fig.  3  which 
shows  speed  versus  time  for  time  up  to   t  =  30  .   The  wave 
speed  was  computed  from  the  relation 


f" 


(1-u)  dx  (38; 


which  is,  however,  only  valid  in  the  steady  state  (compare 
Eq.  (16)).   The  results  in  Fig.  3  confirm  that  a  new  steady 
state,  as  measured  by  changes  in  the  wave  speed,  is  reached 
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with  speed  "^^   8  .   This  agrees  with  Eq ,  (21b)  for   b  =  0.127. 
Curve  2  in  Fig.  3  summarizes  a  different  set  of  results  which 
were  obtained  under  the  same  conditions  as  those  given  by  Curve 
1,  but  where  the  initial  condition  was  given  only  by  Eq .  (22) . 
In  this  case  the  wave  is  seen  to  propagate  with  a  speed  which 
agrees  well  with  the  exact  solution   S  =  5//6      for  the  duration 
of  the  computation. 

The  sensitivity  of  the  results  to  minor  perturbations  in 
the  initial  data  which  was  seen  in  Figs.  2  and  3,  indicates 
that  the  Fisher  Equation  could  be  an  exacting  test  problem  for 
numerical  wave  propagation  studies.   Accordingly,  the  twelve 
numerical  schemes  which  are  listed  in  Table  I  were  compared 
for  this  problem,  as  will  now  be  discussed  with  the  help  of  Figs 

4  and  5.   These  figures  are  also  plots  of  computed  wave  speed 
versus  time.   Here   Ax  =  1   and   At  =  0.2   and  the  initial  data 
was  specified  with  (22)  (steady  propagating  wave  with  speed 

5  =  5//6"  )  .   For  the  schemes  I  to  VI  shown  in  Fig.  4  the  error 
in  the  computed  wave  speed  appears  to  grow,  at  least  initially, 
linearly  with  time.   The  results  of  schemes  VII  to  XII  are 
shown  in  Fig.  5;  error  growth  is  more  modest  so  that  greater 
accuracy  is  associated  with  most  of  these  schemes. 

Notice  that  the  schemes  in  Table  I  also  resemble  several 
other  standard  schemes:   implicit  (Imp),  Crank-Nicolson  (CN) 
and  explicit  (Exp) .   Each  of  these  methods  may  be  formed  from 
Eqs.  (23a, b)  by  replacing   V.   and   W.   with   u!J"^.   in  either 

the  diffusion  or  reaction  term  or  in  both  terms,  and  by  making 
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the  same  choice  of  the  parameters  y,    6,  a,  q   and   6   as  is 
shown  in  Table  I.   The  schemes  in  Eqs.  (23a,  b  and  c)  produce 
the  same  order   At   truncation  errors  as  the  corresponding  Imp, 
CN  and  Exp  schemes  which  are  shown  in  the  Table.   These  error 
terms  are  listed  first  in  the  truncation  error  column  of  the 

Table.   However,  Eqs.  (23)  generates  additional  terms,  propor- 

2 
tional  to   rAt   and   rAt   ,   which  are  given  in  the  Table.   In 

the  case  of  the  standard  explicit  method  Eqs.  (23a,b)  become 

identical  and  reduce  to  the  standard  explicit  scheme. 

The  underestimate  of  the  steady  wave  speed  which  is  seen  in 
Fig.  4  for  the  standard  explicit  method  VI  (a=l,  6=y=q=B=0)    is 
primarily  due  to  a  first  order  numerical  diffusion  error.   This 
is  demonstrated  by  comparing  the  results  of  this  scheme  with 
the  results  of  scheme  XI  (Fig.  5),  which  is  identical  to  scheme 
VI,  except  that  the  physical  diffusion  coefficient  d   was 
increased  by  a  factor  of   1  +  (S  At/2d)  .   Here   S^'^^   was 
computed  explicitly  from  Eq .  (38)  as  the  computation  proceeded. 
This  increase  in   d   minimizes  the  numerical  diffusion  effect 
in  steady  wave  propagation  as  is  seen  in  the  truncation  error 
column  of  scheme  XI  in  Table  I. 

A  penalty  which  is  associated  with  the  significant  improve- 
ment in  the  accuracy  of  the  explicit  method  which  is  seen  in 
Scheme  XI  is  that  the  diffusional  stability  criterion   r  <^  1/2 

(Compare  Eq .  (27)  for   y  =  0  )  becomes  more  restrictive  for 
higher  wave  speeds.   For   S  >>  d/Ax  ,   the  stability  condition 

(27)  is  replaced  by  the  C.F.L.  condition   SAt/Ax  <  1  .   However, 
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scheme  XI  is  second  order  accurate  for  the  case  of  steady  wave 
propagation,  and  it  is  therefore  of  interest  to  compare  its 
performance  with  others  in  Table  I. 

In  scheme  IX  (6=a=l/2,  Y=q=6=0)  the  leading  truncation 

error  is   — ;r—  du   ^  .   In  contrast  to  scheme  VI  this  is  a 
2     xxt 

numerical  dispersion  error  and  not  a  diffusive  error  in  steady 
propagation.   A  comparison  of  scheme  IX  (Fig.  5),  v/ith  the 
results  of  scheme  VI  (Fig.  4)  shows  that  this  first  order  dis- 
persive error  does  not  deteriorate  accuracy  as  much  as  the 
diffusive  error.   In  fact,  the  accuracy  of  the  method  IX  is 
seen  to  be  comparable  to  that  of  the  second  order  scheme  XI . 

The  parameter   y   in  the  diffusion  terms  of  Eqs.  (23a, b) 
is  given  by  Eq .  (26)  in  scheme  XII.   Here   q=6=0   and   6=a=l/2  . 

This  optimal  choice  of  the  parameter   y   results  in  a  leading 

2    y   h 
truncation  error  of  the  form   rAt    L'  ^(u)  .   The  terms 

y    6 
L'   (u)   m  Table  I  are  derived  from  a  Taylor  series  expansion 

of  the  methods  followed  by  an  analysis  similar  to  that  which  led 

to  Eq.  (25).   In  Fig.  5  it  is  seen  that  scheme  XII  performs  as 

well  as  the  second  order  accurate  scheme  XI.   From  this  it  would 

2   Y  6 
be  concluded  that  the   rAt   L''  (u)   truncation  error  terms 

are  small.   The  results  which  correspond  to  schemes  VIII  a,  b 

and  X  are  also  shown  in  Fig.  5,  and  they  are  seen  to  give 

accurate  results.   The  three  schemes  have  the  same  leading 

truncation  error.   However,  differences  in  their  respective 

truncation  errors  do  occur  in  the  higher  order  terms  (see  Table  I] 

and  the  agreement  in  their  results  confirms  that  the  effect  of 


2 
these  0(rAt  )  terms  is  small. 

The  additional  term  ^  (F  +  SF  ) ,  which  vanishes  in  the 
steady  state,  in  scheme  X  (Y=a=l,  6=3=0)  comes  from  evaluating 
the  reaction  rate  explicitly  as   F (U .   .  )   with   q  =  -SAt/2Ax. 

The  reaction  rate  at   x.+qAx   (here,  -1  <  q  <  0)   was  determined 
by  linear  interpolation  between  the  rates  at   x .   and   x .  , 

with   S      given  by  Eq.  (38) .   This  method  of  computing  the 
reaction  rate  minimizes  a  truncation  error,   — j^—  F   ,   which 
arises  from  the  similar  but  less  accurate  scheme  V   (Y=a=l, 
6=q=3=0) .   This  may  be  seen  by  comparing  scheme  V  in  Fig.  4 
with  scheme  X  in  Fig.  5.   The  reaction  rate  is  evaluated  im- 
plicitly in  scheme  Villa  (y=1/  6=a=l/2,  q=3=0)  and  is  computed 
explicitly,  but  with  the  same  leading  truncation  error,  in  scheme 
Vlllb  (y=1,  6=q=0,  a  =  I  ,  3  =  -j)  •   However,  the  methods  Villa 
and  Vlllb  (dashed  curve)  are  seen  to  produce  similar  results 
in  Fig.  5. 

The  methods  VIII  a,  b  and  X  have  a  first  order  numerical 
A  comparison  of  their  results 
with  those  of  the  previously  discussed  scheme  IX,  which  has  a 
truncation  error  — p—  du   .  ,  shows  a  corresponding  change  of 
sign  in  the  error  in  the  computed  wave  speed. 

The  first  order  dispersive  error  is  larger  for  scheme  VII 
and  is  (i^  +  4r)dAtu^^^  (Y  =  2,  q=3=0  and  6=a=i^)  .   Referring  to 
Fig.  5,  this  scheme  is  seen  to  overestimate  the  steady  wave 
speed  significantly.   Therefore,  it  would  be  concluded  that 
largely  numerical  dispersion  errors  can  deteriorate  the  accuracy 
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of  a  method. 

The  results  shown  in  Fig.  4  also  include  schemes  I  to  IV 
which  have  not  yet  been  discussed.   Scheme  I  gave  the  poorest 
results  in  the  study.   This  scheme  has   y  =  2   and  the  reaction 
rate  is  evaluated  fully  implicitly  (6=1,  a=6=q=0) .   The  scheme 
has  a  dispersive  error   4rAtdu   .   which  is  greater  than  that 
of  schemes  Villa,  b,  X,  but  less  than  that  of  scheme  VII. 
However,  the  method  has,  in  addition,  a  numerical  diffusion 
error  ~—  U  .  .   A  comparison  between  schemes  VII,  Villa,  b  and 
X  (Fig.  5)  and  scheme  I  (Fig.  4)  shows  that  this  numerical 
diffusion  term  deteriorates  the  accuracy  of  the  method  greatly. 

Finally,  a  comparison  of  the  performance  of  schemes  I  to 
VI  (Fig.  4)  with  that  of  the  generally  superior  schemes  VII  to 
XII  (Fig.  5)  shows  the  importance  of  an  accurate  evaluation  of 
the  reaction  rate.   Schemes  VII  to  XII  have  first  order  trun- 
cation errors  '\^At  (U  .  -  F.  )  .   This  form  of  truncation  error 
would  lead  to  a  second  order  method  in  the  absence  of  physical 
diffusion,  i.e.  d=0  ,   and  it  is  seen  to  be  preferable  also 
for  the  case   d  7^  0  .   In  the  following  we  present  numerical 
results  which  were  obtained  using  the  explicit  method  Vlllb 
which  satisfies  this  latter  requirement  and  does  not  have  a 
first  order  numerical  diffusion  truncation  error. 

We  next  present  results  for  the  ozone  decomposition  flame 
propagation  problem,  mentioned  earlier,  which  has  also  been 
analysed  in  detail  by  Bledjian  [1] ,  Margolis  [2]  and  by  Hirsch- 
felder  et  al .  [15].   The  elementary  reaction  mechanism  involves 
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atomic  oxygen  0,  molecular  oxygen  O2  and  ozone  0-,  through  the 
reactions 

,  f 


M 


■^   0, 


M 


0+0. 


O2  +  M 


0   +   0  t  M 


In  the  above  reaction  steps  M  stands  for  any  of  the  three 
chemical  species   0,  0^  and  O^  (hence  there  are  7  reactions)  and 
the  superscripts,  f,b,  refer  to  the  forward  and  backward  reactions, 

The  reaction  rates  in  Eqs.  (10)  and  (11)  are  computed  from 


m=l 


in 


(39) 


where 


.   3 

m  n 
n=l 


n,m 


and   R    is  given  correspondingly  with   f   replaced  by   b   and 

with  the  order  constants   a'   replaced  by   a"   and  the  stoichio- 

m  m 


metric  coefficients   v '     replaced  by   v " 

n  ,m  -'         n,m 


These  last  two 


quantities  are  related  through 
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3    -  ..3   ^" 

The  rate  constants  are  given  by  the  Arrhenius  expression 
.     ^   s^   -E^/RT 

with  a  corresponding  expression  for   k   .   The  values  of  the 


Dimension  less  governing  reaction  diffusion  equations  are 
as  given  in  Eqs.  (10)  to  (13)  ,  but  with  p'^V   =  x^/t   , 

X  =  PqX/x^  ,   t  -  t/t^   and  N  =  3.   y^"""^,  y^^^  and  y^^^  are 

the  mass  fractions  of   0  ,  0^      and   0^  ,   and   T,  p,  bi^    '       and 

(k) 
h      are  non-dimensionalized  with   T   ,  p   ,  t    and   p  T  c   , 
■'-  o    o    o        o  o  p 

respectively.   t^  is  a  characteristic  time  scale  chosen  as 

5.878  X  io~   sec,  and  a  length  scale   x   =  5.048  x  10~^  a/cm^ 

is  also  used  in  the  nondimensionalization .   The  computational 

domain  is   0  ^  x  <_  50   and  the  initial  data  is  prescribed  as 

is  shown  in  Table  II.   This  initial  data  and  nondimensionalization 

was  also  used  by  Margolis  [2] . 

Numerical  results  for  this  ozone  flame  problem  are  shown 
in  Fig.  6  using  the  method  Vlllb  of  Table  I.   The  figure  shows 
our  computed  temperature  and  atomic  mass  fraction  profiles  in 
four  separate  graphs,  each  of  which  spans  the  computational 
domain.   The  graphs  show  results  which  correspond  to   t=0  ,  10, 
20  and  30  where  time  increases  from  top  to  bottom  in  the  figure. 
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The  slope  of  the  oblique  line  which  connects  the  reference  point 
in  the  temperature  profile  (T  =  900°k)  of  each  graph  is  propor- 
tional to  the  flame  speed  since  the  steady  speed  is  invariant  under 
the  transformations  (7a,b) .   The  fact  that  the  line  is  straight 
for   t  >_  10   indicates  that  steady  wave  propagation  has  been 
reached.   However,  the  mass  fraction  of  atomic  oxygen  at   x  =  0, 
which  is  diminishing  in  the  successive  frames,  is  still  several 
orders  of  magnitude  larger  than  its  equilibrium  \/alue.   An 
adaptive  grid  was  used  with  only  the  15  computational  points 
which  are  shown  in  each  frame  with   c  =  4   in  Eq.  (34) .   The 
initial  conditions  in  Table  II  are  such  that  the  ignition 
process  was  essentially  confined  to   0.8£x£l.2  .   During 
the  computation  the  minimum  grid  spacing   z  .    changed  from 
.07  at   t  =  0   to  about  .5  for   t  >  4  . 

Details  of  the  flame's  early  transient  development  and 
its  steady  state  structure  are  shown  in  Figs.  7  and  8  respec- 
tively.  The  figures  show  temperature  and  atomic  oxygen,  molecu- 
lar oxygen  and  ozone  mass  fraction  profiles  in  portions  of  the 
computational  domain.   The  solid  curves  show  the  numerical 
results  obtained  by  Margolis  [2],  who  used  an  equally  spaced 
mesh  (Ax  =  0.2)  and  a  fourth  order  accurate  method  of  lines. 
A  comparison  of  the  results  shows  excellent  agreement  of  the 
profiles  in  Fig.  7  at   t  =  0.1  .   In  Fig.  8  (t  =  35)  the  computed 
flame  profiles  also  show  good  agreement,  but  the  flame  location 
differs.   This  is  due  to  a  difference  in  the  computed  steady 
state  wave  speed  of  the  flame  in  the  two  calculations. 
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This  difference  in  flame  speed  can  be  seen  in  Fiq.  9. 
This  figure  shows  how  the  computed  steady  state  flame  speed 
varied  as  a  function  of  the  number  of  grid  points  used  in  five 
separate  computations.   For  the  conditions  of  Fig.  8  (15  grid 
points)  the  computed  wave  speed  differs  from  that  given  by 
Margolis  [2]  (dashed  line  in  Fig.  9)  by  about  2%.   Notice  in 
Fig.  9  that  the  rate  of  convergence  of  the  results  with  an  increase 
in  the  number  of  grid  points  indicates  that  the  present  numerical 
method  appears  to  be  second  order  accurate.   In  these  computa- 
tions the  constant   c   in  Eq .  (34)  was  varied  in  proportion  with 
the  total  number  of  grid  points  as  is  indicated  in  the  figure. 

The  flame  speed  was  computed  using  the  method  which  was 
indicated  in  the  discussion  of  Fig.  6.   It  was  found  to  vary 
by  less  than  1%  in  each  of  the  11,  15,  21  and  30  grid  point  cal- 
culations for   t  ^  10.   The  average  speed  is  given  in  Fig.  9 
for   10  £  t  <_  35   and  has  the  value   S  =  49.8  ±  0.1  cm/s   for 
30  grid  points. 

The  computer  execution  time  (CDC  6600)  is  also  shown  in 
Fig.  9.   It  increases  nearly  linearly  with  the  number  of  grid 
points.   No  program  optimization  was  performed.   Each  flame 
calculation  required  about  7,500  computational  cycles  to  reach 
t  =  35.   The  time  step  ranged  from   At  =  0.0025    in  the  transient 
to  an  average  value  of   At  %  0.005  ,   or   At  '^^  3ys   for   t  >  4. 
Time  step  values  were  not  given  by  Bledjian  [1]  or  Margolis  [2] , 
but  the  values  given  by  Lund  [4],  for  a  similar  ozone  model 
problem,  are  comparable;   At  '^^  0.007,   or   At  ^i   4ijs.   In  his 
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implicit  code  the  timestep  was  controlled  by  specifying  a  maximum 
allowable  number  of  iterations  for  convergence  of  the  solution  in 
each  computational  cycle. 

The  agreement  between  the  timestep  of  the  present  explicit 
method  with  that  of  the  implicit  method  of  Lund  [4]  could  imply 
that  each  of  the  mass  fraction  and  temperature  reaction-diffusion 
equations  has  a  similar  degree  of  stiffness  for  this  ozone  flame 
problem.   In  other  flame  problems,  component  species  may  be  present 
which  do  not  influence  the  eneray  release  and/or  flame  speed 
significantly,  but  v;hich  do  limit  the  timestep.   In  this  case 
the  use  of  stiffly  stable  methods  for  the  reaction  rate  would  be 
advantageous.   An  explicit,  stiffly  stable  method  for  the  reaction 
rate  is  currently  being  developed  and  tested  for  this  case.   In 
the  method  the  reaction  terms  in  Eqs.  (23a,  b)  are  replaced  by 

,n+l 


^i   K,r 


^'K,:> 


<,1 


"ko 


where  the  superscripts   f,b   refer  to  the  positive  and  negative 
contributions  to  the  overall  rate   F  =  F  -F    (Compare  Eq .  (39)). 

The  effect  of  this  method  is  illustrated  by  applying  it 
to  the  Fisher  Equation  (37)  with  d  =  0  which  has  the  exact 
solution  / 

u(t)  =   1/  1  +  (-^U"^ 


where   u    is   u(t  =  0) .   The  numerical  solution  to  the  scheme 
o 


At 
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can   be    shown    to   be 


/I'  •  ^ 


(1    +    At) 
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5.  CONCLUSIONS 

The  model  equation  parameter  study  in  section  4  has  shown 
that  numerical  diffusion  effects  can  reduce  the  accuracy 
of  a  numerical  method  significantly.   Numerical  dispersion 
truncation  errors  were  also  found  to  be  detrimental  to  accuracy, 
but  only  if  they  were  sufficiently  large.   The  study  also  points 
out  that  the  reaction  rate  term  should  be  accurately  represented 
in  the  finite  difference  scheme.   Generally,  methods  which  were 
second  order  accurate  for  the  case  of  no  physical  diffusion, 
d  =  0,   were  found  to  be  best. 

The  explicit  scheme  Vlllb,  which  was  selected  for  the  ozone 
flame  calculations,  was  found  to  behave  like  a  second  order 
accurate  scheme.   The  scheme  has  no  first  order  numerical  dif- 
fusion error  and  from  the  results  it  is  concluded  that  the 

3   2 
effect  of  the  dispersive  error  term  and  the   0(At  /h  )   trunca- 
tion terms  was  minimal. 

The  favorable  comparison  of  the  ozone  flame  results  with 
the  results  of  the  fourth  order  accurate  method  of  Margolis  [2] 
indicates  that  a  second  order  method  is  sufficiently  accurate 
for  practical  computations.   Finally,  the  fact  that  the  numerical 
timestep  of  the  present  explicit  method  was  found  to  be  comparable 
to  the  timestep  used  in  the  implicit  method  of  Lund  [4]  indicates 
that  explicit  numerical  methods  could  prove  to  be  competitive 
in  certain  combustion  problems. 
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FIGURE  AND  TABLE  CAPTIONS 


FIGURE  1  -  Schematic  Diagram  of  One-dimensional  Flame  Propagation 
Problem. 

FIGURE  2  -  Numerical  Solution  of  Fisher's  Equation  with  Initial 
Data- Equation  (22)  (U  >  5  x  10"^)  and  Equation  (21a) 
(U  <_  5  X  10-6)  ^ 

FIGURE  3  -  Computed  Wave  Speed  versus  Time  for  Fisher's  Equation. 
Initial  Data  -  Curve  1,  Equations  (22)  and  (21a)  - 
Curve  2,  Equation  (22). 

FIGURE  4  -  Computed  Wave  Speed  versus  Time  for  Fisher's  Equation. 
Initial  Data  -  Equation  (22).   Schemes  I  -  VI. 

FIGURE  5  -  Computed  Wave  Speed  versus  Time  for  Fisher's  Equation. 
Initial  Data  -  Equation  (22),   Schemes  VII-XII. 

FIGURE  6  -  Profiles  of  Temperature  D  and  Atomic  Oxygen  Mass 

Fraction   o   in  the  Propagatina  Ozone  Decomposition 
Flame  for   t  =  0,  10,  20  and  30.   15  grid  points, 
c  =  4   in  Equation  (34). 

FIGURE  7  -  Profiles  of  Temperature  n  and  Atomic  Oxygen   o  , 
Molecular  Oxygen  X    ,    Ozone   A   mass  fractions  for 
t  =  0.1  .   15  grid  points,   c  =  4   in  Equation  (34) . 

—  Results  of  Margolis  [2]  . 

FIGURE  8  -  Profiles  of  Temperature  n   and  Atomic  Oxygen  o  , 

Molecular  Oxygen  X  i      Ozone   A   mass  fractions  for 
t  =  35  ,   15  grid  points,   c  =  4   in  Equation  (34) . 

—  Results  of  Margolis  [2]. 

FIGURE  9  -  Variation  of  Computed  Steady  Wave  Speed  with  Number 
of  Grid  Points  M.   c-constant  in  Equation  (34). 

—  Wave  Speed  of  Margolis  [2] . 

TABLE  I   -   Parameter  Study  of  Numerical  Method,  Equations  (23a, 
b  and  c) .   (Imp,  CN,  Exp  -  Analogous  Implicit, 
Crank-Nicolson  and  Explicit  Methods) . 

TABLE  II  -  Constants  and  Initial  Data  for  Ozone  Decomposition 
Flame  Studies. 
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TABLE    II 


Rate  Data 

m=l,2,3 

m=4 

m=5,6,7 

E  (cal/mole) 
m 

2.414  X  10^ 

6.000  X  10^ 

1.1735  X  10^ 

B^(l/sec°K  "^) 

6.760  X  10^ 

4.580  X  10^ 

5.710  X  10^ 

<<-) 

2.5 

2.5 

2.5 

E^ (cal/mole) 

0 

9.921  X  10^ 

0 

B^l/sec°K  "^) 
m 

1.180  X  10^ 

1.880  X  10^ 

2.470  X  10^ 

=^-) 

3.5 

2.5 

3.5 

Thermodynamic  Data 

k=l 

k=2 

k  =  3 

h^"^^  (cal/mole) 

5.8675  X  lo"^ 

0 

3.4535  X  10^ 

W^^^  (g/mole) 

16 

32 

48 

Cp(cal/g°K) 

0.2524 

p  (atm.) 

0.821 

T^(OK) 

300 

p^(a/cm  ) 

1.201  X  10"^ 

t^(sec) 

5.878  X  lo"^ 

x^(g/cm  ) 

5.048  X  lo"^ 

Initial  Data 
g(x)  =  cos^[l[^)'' 

0  <  x  < 

1.2 

1.2 

<  X  <_  5  0 

Y^^^  (x,0  ) 

5  X  10"^g(x) 

0 

Y^2)  (x,0) 

(2+g(x))/3-Y 

^^^x,0) 

2/3 

y(^)(x,0) 

1-Y^^^  (x,0)- 

y(2)(x,o) 

1/3 

T(x,0) 

1  +  19g(x 

)/6 
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